Project Summary

This project aimed to look at the spatial variability of color morph and Symbiodinium clades C and D present in the Kane’ohe Bay, O’ahu, Hawai’i population of Montipora capitata. We investigated the distributions at scales ranging from location within the bay to location on an individual reef. We also looked at differences among reef types (fringing vs. patch) and depths. Heterogeneous mixtures of symbiont clades were considered in the analysis for spatial patterns. By investigating spatial variability of Symbiodinium, we furthered the understanding of potential stress-response in Kane’ohe Bay.

Library Packages

knitr::opts_knit$set(root.dir = normalizePath(".."))
# setwd("~/Symcap") #This will not work on other people's computers and you should not need this - Rproj automatically sets wd to root of repository
# I don't thing all of these packages are actually being used in the code...
library(data.table)
library(lsmeans)
library(devtools)
library(plyr)
library(reshape2)
library(popbio)
library(RgoogleMaps)
library(plotrix)
library(zoo)
library(rgdal)
library(car)
library(scales)
library(png)
library(pixmap)
library(ecodist)
library(cluster)
library(fpc)
library(clustsig)
library(mapplots)
library(foreign)
library(nnet)
library(ggplot2)
library(mlogit)

Import/Merge Collection and qPCR Data

# Import sample metadata
Coral_Data <- read.csv("Data/Collection Data/Coral_Collection.csv")
Coral_Data$Depth..m. <- as.numeric(as.character(Coral_Data$Depth..m.))

# Import qPCR data
source_url("https://raw.githubusercontent.com/jrcunning/steponeR/master/steponeR.R")
Mcap.plates <- list.files(path = "Data/qPCR_data", pattern = "txt$", full.names = T)
Mcap <- steponeR(files = Mcap.plates, delim="\t",
                 target.ratios=c("C.D"), 
                 fluor.norm = list(C=2.26827, D=0), 
                 copy.number=list(C=33, D=3))
Mcap <- Mcap$result
## Remove positive and negative controls from qPCR data
Mcap <- Mcap[grep("+", Mcap$Sample.Name, fixed=T, invert = T), ]
Mcap <- Mcap[grep("NTC", Mcap$Sample.Name, fixed = T, invert = T), ]
Mcap <- Mcap[grep("PCT", Mcap$Sample.Name, fixed = T, invert = T), ]
## Change "Sample.Name" column to "Colony"
colnames(Mcap)[which(colnames(Mcap)=="Sample.Name")] <- "Colony"
## Filter out detections without two technical replicates
Mcap$fail <- ifelse(Mcap$C.reps < 2 & Mcap$D.reps < 2, TRUE, FALSE)
fails <- Mcap[Mcap$fail==TRUE, ]
Mcap <- Mcap[which(Mcap$fail==FALSE),]
## Set C:D ratio to -Inf (when only D) and Inf (when only C)
Mcap$C.D[which(Mcap$C.reps<2)] <- -Inf
Mcap$C.D[which(Mcap$D.reps<2)] <- Inf
## Order qPCR data by Colony
Mcap <- Mcap[with(Mcap, order(Colony)), ]
## Calculate proportion C and proportion D
Mcap$propC <- Mcap$C.D / (Mcap$C.D+1)
Mcap$propD <- 1-Mcap$propC
Mcap$propD[which(Mcap$C.D==-Inf)] <-1
Mcap$propC[which(Mcap$C.D==-Inf)] <-0
Mcap$propD[which(Mcap$C.D==Inf)] <-0
Mcap$propC[which(Mcap$C.D==Inf)] <-1
## Categorize each sample as C-dominated or D-dominated
Mcap$Dom <- ifelse(Mcap$propC>Mcap$propD, "C", "D")

# Merge sample metadata with qPCR data
Symcap<-merge(Coral_Data, Mcap, by="Colony", all=T)
Symcap <- Symcap[with(Symcap, order(Colony)), ]
## Define category for symbiont mixture
Symcap$Mix <- factor(ifelse(Symcap$propC>Symcap$propD, ifelse(Symcap$propD!=0, "CD", "C"), ifelse(Symcap$propD>Symcap$propC, ifelse(Symcap$propC!=0, "DC", "D"), NA)), levels = c("C", "CD", "DC", "D"))
## Define reef.area category as either Top or Slope
Symcap$Reef.Area <- ifelse(Symcap$Reef.Area!="Top", yes = "Slope", no = "Top")

# Import color data from each observer
colscores <- read.csv("Data/Collection Data/Color_Scores.csv")
## Determine color scores by majority rule, and number of agreements
colscores$majority <- apply(colscores[,-1], 1, function(x) names(table(x))[which.max(table(x))])
colscores$n.agree <- apply(colscores[,-c(1,7)], 1, function(x) max(table(x)))
## Set which set of color scores will be used in analysis (=majority)
Symcap <- merge(Symcap, colscores[,c("Colony","majority","n.agree")], all=T)
Symcap$Color.Morph <- Symcap$majority

Change Bay Area for 18, 26 and F7-18

# PROVIDE RATIONALE FOR MAKING THESE CHANGES
Symcap$Bay.Area[which(Symcap$Reef.ID=="26")] <- "Central"
Symcap$Bay.Area[which(Symcap$Reef.ID=="18")] <- "Southern"
Symcap$Bay.Area[which(Symcap$Reef.ID=="F7-18")] <- "Southern"

Adjust Depth by Mean Sea Level

To account for the influence of tides, depth was adjusted according to the differences in mean sea level from the recorded sea level on each collection day to the nearest 6-minute interval. Mean sea level was obtained from NOAA daily tide tables for Moku o Lo’e.

JuneTide=read.csv("Tide Tables/Station_1612480_tide_ht_20160601-20160630.csv")
JulyTide=read.csv("Tide Tables/Station_1612480_tide_ht_20160701-20160731.csv")
AugustTide=read.csv("Tide Tables/Station_1612480_tide_ht_20160801-20160812.csv")
Tide<-rbind(JuneTide, JulyTide, AugustTide)
Tide$Time <- as.POSIXct(Tide$TimeUTC, format="%Y-%m-%d %H:%M:%S", tz="UTC")
attributes(Tide$Time)$tzone <- "Pacific/Honolulu"
Symcap$Time2 <- as.POSIXct(paste(as.character(Symcap$Date), as.character(Symcap$Time)),
                                format="%m/%d/%y %H:%M", tz="Pacific/Honolulu")
Symcap$Time=Symcap$Time2

# Add estimated times for missing values
Symcap$Time[170] <- as.POSIXct("2016-06-14 12:07:00")
Symcap$Time[177] <- as.POSIXct("2016-06-14 12:20:00")
Symcap$Time[178] <- as.POSIXct("2016-06-14 12:22:00")
Symcap$Time[180] <- as.POSIXct("2016-06-14 13:08:00")
Symcap$Time[187] <- as.POSIXct("2016-06-14 12:42:00")
Symcap$Time[188] <- as.POSIXct("2016-06-14 12:44:00")
Symcap$Time[206] <- as.POSIXct("2016-06-16 13:10:00")
Symcap$Time[208] <- as.POSIXct("2016-06-16 13:24:00")
Symcap$Time[211] <- as.POSIXct("2016-06-16 12:37:00")
Symcap$Time[218] <- as.POSIXct("2016-06-16 12:27:00")
Symcap$Time[448] <- as.POSIXct("2016-07-16 13:32:00")

Round6 <- function (times)  {
  x <- as.POSIXlt(times)
  mins <- x$min
  mult <- mins %/% 6
  remain <- mins %% 6
  if(remain > 3L) {
    mult <- mult + 1
  } else {
    x$min <- 6 * mult
  }
  x <- as.POSIXct(x)
  x
}

Symcap$Time.r <- Round6(Symcap$Time)
Tide$Time.r <- Tide$Time
merged<-merge(Symcap, Tide, by="Time.r", all.x=T)
merged$newDepth <- merged$Depth..m.- merged$TideHT

Coral Collection Map

Collection sites represented on this map indicate the 9 fringe and 16 patch reef sites for a total of 25 collection reefs. In total, 707 colonies were collected for the analysis.

# Define map area and get satellite data from Google using RgoogleMaps
KB <- c(21.46087401, -157.809907) 
KBMap <- GetMap(center = KB, zoom = 13, maptype = "satellite", SCALE = 2, GRAYSCALE = FALSE)
# WHAT IS THE PURPOSE OF GENERATING MAP, SAVING TO FILE, THEN LOADING FROM FILE?
save(KBMap, file = "KBMap.Rdata")
load("KBMap.Rdata") 

# Get representative coordinates for each reef based on sample GPS data
Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
newcoords <- LatLon2XY.centered(KBMap, XY$Latitude, XY$Longitude, zoom=13)
XY$X <- newcoords$newX
XY$Y <- newcoords$newY
XY <- subset(XY, Reef.ID!="37")  # WHAT IS THIS DOING? WHY?

# Plot KBay Map
par(oma=c(3,3,0,0))
PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude, col=153, pch=21, bg="#FF7F50", lwd=2, cex = 1.2)
axis(1, at = LatLon2XY.centered(KBMap, NA, c(-157.85, -157.81, -157.77))$newX, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("157.85°W", "157.81°W", "157.77°W"), padj = -2.5, cex.axis = 0.75)
axis(2, at = LatLon2XY.centered(KBMap, c(21.42, 21.46, 21.50), NA)$newY, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("21.42°N", "21.46°N", "21.50°N"), padj = 0.5, hadj = 0.60, las = 1, cex.axis = 0.75)
par(new=T, mar=c(14,21,0,0))
HI <- readOGR("coast_n83.shp", "coast_n83") 
HI <- spTransform(HI, CRS("+proj=longlat +datum=NAD83")) 
plot(HI, xlim=c(-158.3, -157.6), ylim=c(21.35, 21.6), lwd=0.4, col="gray", bg="white")
rect(-157.87, 21.41, -157.75, 21.52)
box()

Symbiont Dominance and Color Morph Across Geographic Scales

Significant Differences

Dominant Symbiont per Reef Area

## I DON'T SEE ANYWHERE IN THE MANUSCRIPT THAT ANY EFFECT OF "REEF.AREA" (=TOP VS SLOPE) IS REPORTED OR DISCUSSED. IF NOT IN MANUSCRIPT, DELETE THIS ANALYSIS FROM MARKDOWN.
results <- table(Symcap$Dom, Symcap$Reef.Area)
chisq.test(results)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  results
## X-squared = 136.26, df = 1, p-value < 2.2e-16
prop.table(results, margin = 2)
##    
##         Slope       Top
##   C 0.7767857 0.3294574
##   D 0.2232143 0.6705426
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray100"), xlab = "Reef Area", ylab = "Proportion of colonies")
legend("topright", legend=c("C-dom.", "D-dom."), fill=c("gray10", "gray100"), inset = c(-.2, 0), xpd = NA)

Dominant Symbiont per Color Morph

results=table(Symcap$Dom, Symcap$Color.Morph)
chisq.test(results)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  results
## X-squared = 263.61, df = 1, p-value < 2.2e-16
prop.table(results, margin = 2)
##    
##          Brown     Orange
##   C 0.92352941 0.32513661
##   D 0.07647059 0.67486339
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray100"), xlab = "Color Morph", ylab = "Symbiont Proportion")
legend("topright", legend=c("C", "D"), fill=c("gray10", "gray100"), inset = c(-.2, 0), xpd = NA)

Symbiont Community Composition per Color Morph

results=table(Symcap$Mix, Symcap$Color.Morph)
chisq.test(results)
## 
##  Pearson's Chi-squared test
## 
## data:  results
## X-squared = 266.44, df = 3, p-value < 2.2e-16
prop.table(results, margin = 2)
##     
##            Brown      Orange
##   C  0.794117647 0.286885246
##   CD 0.129411765 0.038251366
##   DC 0.073529412 0.653005464
##   D  0.002941176 0.021857923
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray50", "gray85", "gray100"), xlab = "Color Morph", ylab = "Symbiont Community Composition")
legend("topright", legend=c("C", "CD", "DC", "D"), fill=c("gray10", "gray50", "gray85", "gray100"), inset = c(-.2, 0), xpd = NA)

Symbiont Community Composition per Reef Area

Symcap$Reef.Area <- ifelse(Symcap$Reef.Area!="Top", yes = "Slope", no = "Top")
results=table(Symcap$Mix, Symcap$Reef.Area)
chisq.test(results)
## 
##  Pearson's Chi-squared test
## 
## data:  results
## X-squared = 138.97, df = 3, p-value < 2.2e-16
prop.table(results, margin = 2)
##     
##            Slope         Top
##   C  0.678571429 0.275193798
##   CD 0.098214286 0.054263566
##   DC 0.214285714 0.651162791
##   D  0.008928571 0.019379845
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray50", "gray85", "gray100"), xlab = "Reef Area", ylab = "Symbiont Community Composition")
legend("topright", legend=c("C", "CD", "DC", "D"), fill=c("gray10", "gray50", "gray85", "gray100"), inset = c(-.2, 0), xpd = NA)

Symbiont Community Composition per Dominant Symbiont

results=table(Symcap$Mix, Symcap$Dom)
chisq.test(results)
## 
##  Pearson's Chi-squared test
## 
## data:  results
## X-squared = 707, df = 3, p-value < 2.2e-16
prop.table(results, margin = 2)
##     
##               C          D
##   C  0.86635945 0.00000000
##   CD 0.13364055 0.00000000
##   DC 0.00000000 0.96703297
##   D  0.00000000 0.03296703
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray 10", "gray 85", "gray 40", "gray100"), xlab = "Dominant Symbiont", ylab = "Symbiont Mixture Proportion")
legend("topright", legend=c("C", "CD", "DC", "D"), fill=c("gray10", "gray85", "gray40", "gray100"), inset = c(-.2, 0), xpd = NA)

When looking at D-only colonies, the Ct values are on the higher end (35 or greater) indicating poor amplification or the potential for C amplification being pushed back in the cycle. This is supported by the fact that 5 of the 9 D-only colonies had C present in 1 qPCR replicate.

Proportion of D When Present in Mixture

propD <- merged$propD[which(merged$propD > 0 & merged$propD < 1)]
hist(propD, xlab = "Proportion of Clade D", ylab = "Number of Samples", main = "", col = "gray75")

Color Morph per Reef Area

Symcap$Reef.Area <- ifelse(Symcap$Reef.Area!="Top", yes = "Slope", no = "Top")
results=table(Symcap$Color.Morph, Symcap$Reef.Area)
chisq.test(results)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  results
## X-squared = 57.683, df = 1, p-value = 3.08e-14
prop.table(results, margin = 2)
##         
##              Slope       Top
##   Brown  0.5902004 0.2906977
##   Orange 0.4097996 0.7093023
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c("gray10", "gray100"), xlab = "Reef Area", ylab = "Color Morph Proportion")
legend("topright", legend=c("brown", "orange"), fill=c("gray10", "gray100"), inset = c(-.2, 0), xpd = NA)

Color Morph per Bay Area

results=table(merged$Color.Morph, merged$Bay.Area)
prop.table(results, margin = 2)
##         
##            Central  Northern  Southern
##   Brown  0.5288462 0.4358974 0.4769737
##   Orange 0.4711538 0.5641026 0.5230263
chisq.test(results)
## 
##  Pearson's Chi-squared test
## 
## data:  results
## X-squared = 3.5162, df = 2, p-value = 0.1724

Color Morph per Reef ID

results=table(Symcap$Color.Morph, Symcap$Reef.ID)
chisq.test(results)
## 
##  Pearson's Chi-squared test
## 
## data:  results
## X-squared = 34.081, df = 24, p-value = 0.08324
Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
propCol=table(Symcap$Color.Morph, Symcap$Reef.ID)
propCol=prop.table(propCol, margin = 2)
propCol <- as.data.frame.matrix(propCol)
props <- data.frame(t(propCol))
props$Reef.ID <- rownames(props)
XY<-merge(XY, props, by="Reef.ID", all=T)
newcoords <- LatLon2XY.centered(KBMap, XY$Latitude, XY$Longitude, zoom=13)
XY$X <- newcoords$newX
XY$Y <- newcoords$newY
PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude)

rownames(XY) <- XY$Reef.ID
XY <- XY[, -1]
XY <- na.omit(XY)
apply(XY, MARGIN=1, FUN=function(reef) {
  floating.pie(xpos = reef["X"], ypos = reef["Y"], 
               x=c(reef["Orange"], reef["Brown"]), radius = 7, col = c("orange", "sienna"))
})

Mantel Test for Spatial Autocorrelation

Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
propCol=table(Symcap$Color.Morph, Symcap$Reef.ID)
propCol=prop.table(propCol, margin = 2)
propCol <- as.data.frame.matrix(propCol)
props <- data.frame(t(propCol))
props$Reef.ID <- rownames(props)
XY<-merge(XY, props, by="Reef.ID", all=T)
reef.dists <- dist(cbind(XY$Longitude, XY$Latitude))
col.dists <- dist(XY$Brown)
set.seed(12456)
mantel(col.dists~reef.dists)
##     mantelr       pval1       pval2       pval3   llim.2.5%  ulim.97.5% 
## -0.04784488  0.66800000  0.33300000  0.59100000 -0.15930224  0.01286751

Color Morph per Bay Area Adjusted for Depth

To discount the vertical spatial plane (depth), these results are adjusted to “assume” that all colonies were collected from the same depth, therefore considering simply the coordinate plane as the distributional influence on the resulting values.

merged$Color <- ifelse(merged$Color.Morph=="Orange", 1, 0)
dat <- aggregate(data.frame(prop=merged$Color), by=list(Reef.ID=merged$Reef.ID), FUN=mean, na.rm=T) 
mod <- glm(Color ~ newDepth + Reef.ID, data=merged)
library(lsmeans)
lsm <- lsmeans(mod, specs="Reef.ID")
res <- merge(dat, data.frame(summary(lsm))[,c(1:2)], by="Reef.ID")
colnames(res) <- c("Reef.ID", "Raw", "O.Adj")

Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
propCol=table(Symcap$Color.Morph, Symcap$Reef.ID)
propCol=prop.table(propCol, margin = 2)
propCol <- as.data.frame.matrix(propCol)
props <- data.frame(t(propCol))
props$Reef.ID <- rownames(props)
XY<-merge(XY, props, by="Reef.ID", all=T)
newcoords <- LatLon2XY.centered(KBMap, XY$Latitude, XY$Longitude, zoom=13)
XY$X <- newcoords$newX
XY$Y <- newcoords$newY
PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude)
XY2 <- merge(XY, res, by = "Reef.ID")
XY2$B.Adj <- 1-XY2$O.Adj

PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude)

rownames(XY2) <- XY2$Reef.ID
XY2 <- XY2[, -1]
XY2 <- na.omit(XY2)
apply(XY2, MARGIN=1, FUN=function(reef) {
  floating.pie(xpos = reef["X"], ypos = reef["Y"], 
               x=c(reef["O.Adj"], reef["B.Adj"]), radius = 7, 
               col = c("orange", "sienna"))
})

Mantel Test for Spatial Autocorrelation

reef.dists <- dist(cbind(XY$Longitude, XY$Latitude))
col.dists <- dist(XY2$O.Adj)
set.seed(12456)
mantel(col.dists~reef.dists)
##     mantelr       pval1       pval2       pval3   llim.2.5%  ulim.97.5% 
##  0.02284001  0.33100000  0.67000000  0.76900000 -0.05447766  0.09818949

Symbiont Community Composition per Bay Area

results=table(Symcap$Mix, Symcap$Bay.Area)
prop.table(results, margin = 2)
##     
##         Central   Northern   Southern
##   C  0.50961538 0.61538462 0.49174917
##   CD 0.08173077 0.05641026 0.09900990
##   DC 0.39423077 0.30256410 0.40594059
##   D  0.01442308 0.02564103 0.00330033
chisq.test(results)
## 
##  Pearson's Chi-squared test
## 
## data:  results
## X-squared = 14.719, df = 6, p-value = 0.02256

Non-Significant Differences

Dominant Symbiont per Reef Type

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  results
## X-squared = 0.77593, df = 1, p-value = 0.3784

Dominant Symbiont per Bay Area

results=table(Symcap$Dom, Symcap$Bay.Area)
chisq.test(results)
## 
##  Pearson's Chi-squared test
## 
## data:  results
## X-squared = 3.8852, df = 2, p-value = 0.1433

Dominant Symbiont per Reef ID

results=table(Symcap$Dom, Symcap$Reef.ID)
chisq.test(results)
## 
##  Pearson's Chi-squared test
## 
## data:  results
## X-squared = 34.598, df = 24, p-value = 0.07459
Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
propDom=table(Symcap$Dom, Symcap$Reef.ID)
propDom=prop.table(propDom, margin = 2)
propDom <- as.data.frame.matrix(propDom)
props <- data.frame(t(propDom))
props$Reef.ID <- rownames(props)
XY<-merge(XY, props, by="Reef.ID", all=T)
newcoords <- LatLon2XY.centered(KBMap, XY$Latitude, XY$Longitude, zoom=13)
XY$X <- newcoords$newX
XY$Y <- newcoords$newY
PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude)

rownames(XY) <- XY$Reef.ID
XY <- XY[, -1]
XY <- na.omit(XY)
apply(XY, MARGIN=1, FUN=function(reef) {
  floating.pie(xpos = reef["X"], ypos = reef["Y"], 
               x=c(reef["C"], reef["D"]), radius = 7, col = c("blue", "red"))
})

Mantel Test for Spatial Autocorrelation

Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
propDom=table(Symcap$Dom, Symcap$Reef.ID)
propDom=prop.table(propDom, margin = 2)
propDom <- as.data.frame.matrix(propDom)
props <- data.frame(t(propDom))
props$Reef.ID <- rownames(props)
XY<-merge(XY, props, by="Reef.ID", all=T)
reef.dists <- dist(cbind(XY$Longitude, XY$Latitude))
dom.dists <- dist(XY$C)
set.seed(12456)
mantel(dom.dists~reef.dists)
##    mantelr      pval1      pval2      pval3  llim.2.5% ulim.97.5% 
## 0.15974663 0.04000000 0.96100000 0.04500000 0.08875833 0.24328773

pval1 indicates that there is a positive correlation between distance and difference among reefs in terms of dominant symbiont clade. Reefs that are closer in proximity to each other are more similar in terms of symbiont dominance. A positive mantelr value indicates a positive correlation.

Clustering of Similar Reefs

Though insignificant differences resulted when investigating symbiont dominance across reefs, reefs closer in proximity appear more similar than reefs further apart. Clustering was tested here.

doms <- hclust(dom.dists)
plot(doms, labels = XY$Reef.ID)

a <- hclust(reef.dists*dom.dists)
plot(a, labels = XY$Reef.ID)

Dominant Symbiont Plots

Raw Values
Latitude=aggregate(Latitude~Reef.ID, data=Symcap, FUN = mean)
Longitude=aggregate(Longitude~Reef.ID, data = Symcap, FUN=mean)
XY<-merge(Latitude, Longitude, by="Reef.ID", all=T)
newcoords <- LatLon2XY.centered(KBMap, XY$Latitude, XY$Longitude, zoom=13)
XY$X <- newcoords$newX
XY$Y <- newcoords$newY
XY <- subset(XY, Reef.ID!="37")

merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
dat <- aggregate(data.frame(prop=merged$Dominant), by=list(Reef.ID=merged$Reef.ID), FUN=mean, na.rm=T) #prop D/reef
mod <- glm(Dominant ~ newDepth + Reef.ID, data=merged)
mod2 <- glm(Dominant ~ newDepth * Reef.ID, data=merged)
library(lsmeans)
lsm <- lsmeans(mod, specs="Reef.ID")
lsm2 <- lsmeans(mod2, specs="Reef.ID")
res <- merge(dat, data.frame(summary(lsm))[,c(1:2)], by="Reef.ID")
res <- merge(res, data.frame(summary(lsm2))[,c(1:2)], by="Reef.ID")
colnames(res) <- c("Reef.ID", "raw.c", "c.adj", "c.adj.int")

XY2 <- merge(XY, res, by = "Reef.ID")
XY2$d.adj <- 1-XY2$c.adj
XY2$d.adj.int <- 1-XY2$c.adj.int
XY2$raw.d <- 1-XY2$raw.c

PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude)

rownames(XY2) <- XY2$Reef.ID
XY2 <- XY2[, -1]
XY2 <- na.omit(XY2)
apply(XY2, MARGIN=1, FUN=function(reef) {
  floating.pie(xpos = reef["X"], ypos = reef["Y"], 
               x=c(reef["raw.c"], reef["raw.d"]), radius = 7, 
               col = c("red", "blue"))
})

Depth Interaction Adjustment
merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
dat <- aggregate(data.frame(prop=merged$Dominant), by=list(Reef.ID=merged$Reef.ID), FUN=mean, na.rm=T) #prop D/reef
mod <- glm(Dominant ~ newDepth + Reef.ID, data=merged)
mod2 <- glm(Dominant ~ newDepth * Reef.ID, data=merged)
library(lsmeans)
lsm <- lsmeans(mod, specs="Reef.ID")
lsm2 <- lsmeans(mod2, specs="Reef.ID")
res <- merge(dat, data.frame(summary(lsm))[,c(1:2)], by="Reef.ID")
res <- merge(res, data.frame(summary(lsm2))[,c(1:2)], by="Reef.ID")
colnames(res) <- c("Reef.ID", "raw.c", "c.adj", "c.adj.int")

XY2 <- merge(XY, res, by = "Reef.ID")
XY2$d.adj <- 1-XY2$c.adj
XY2$d.adj.int <- 1-XY2$c.adj.int
XY2$raw.d <- 1-XY2$raw.c

PlotOnStaticMap(KBMap, XY$Latitude, XY$Longitude)

rownames(XY2) <- XY2$Reef.ID
XY2 <- XY2[, -1]
XY2 <- na.omit(XY2)
apply(XY2, MARGIN=1, FUN=function(reef) {
  floating.pie(xpos = reef["X"], ypos = reef["Y"], 
               x=c(reef["c.adj.int"], reef["d.adj.int"]), radius = 7, 
               col = c("red", "blue"))
})

Color Morph and Dominant Symbiont Interaction Plots

A multinomial logistic regression was performed on the interaction of color morph and dominant symbiont to discount the influence of depth on spatial distribution.

merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
dat <- aggregate(data.frame(prop=merged$Dominant), by=list(Reef.ID=merged$Reef.ID), 
                 FUN=mean, na.rm=T) #prop D/reef
mod <- glm(Dominant ~ newDepth + Reef.ID, data=merged)
mod2 <- glm(Dominant ~ newDepth * Reef.ID, data=merged)
library(lsmeans)
lsm <- lsmeans(mod, specs="Reef.ID")
lsm2 <- lsmeans(mod2, specs="Reef.ID")
res <- merge(dat, data.frame(summary(lsm))[,c(1:2)], by="Reef.ID")
res <- merge(res, data.frame(summary(lsm2))[,c(1:2)], by="Reef.ID")
colnames(res) <- c("Reef.ID", "raw.c", "c.adj", "c.adj.int")

XY2 <- merge(XY, res, by = "Reef.ID")
XY2$d.adj <- 1-XY2$c.adj
XY2$d.adj.int <- 1-XY2$c.adj.int
XY2$raw.d <- 1-XY2$raw.c

manteltable = table(merged$Dom, merged$Color.Morph, merged$Reef.ID)
nc <- aggregate(interaction(merged$Color.Morph, merged$Dom), 
                by=list(merged$Reef.ID), FUN=table)
nc <- data.frame(Reef.ID=as.character(nc[,1]), prop.table(nc[,2], margin=1))
XY3 <- merge(XY2, nc, by="Reef.ID", all=T)

merged$ColDom <- interaction(merged$Color.Morph, merged$Dom)
model <- multinom(ColDom~newDepth+Reef.ID, merged)
means <- lsmeans(model, specs = "Reef.ID")
means
## Warning in sqrt(.qf.non0(object@V, x)): NaNs produced

## Warning in sqrt(.qf.non0(object@V, x)): NaNs produced

## Warning in sqrt(.qf.non0(object@V, x)): NaNs produced

## Warning in sqrt(.qf.non0(object@V, x)): NaNs produced

## Warning in sqrt(.qf.non0(object@V, x)): NaNs produced

## Warning in sqrt(.qf.non0(object@V, x)): NaNs produced

## Warning in sqrt(.qf.non0(object@V, x)): NaNs produced

## Warning in sqrt(.qf.non0(object@V, x)): NaNs produced

## Warning in sqrt(.qf.non0(object@V, x)): NaNs produced

## Warning in sqrt(.qf.non0(object@V, x)): NaNs produced

## Warning in sqrt(.qf.non0(object@V, x)): NaNs produced

## Warning in sqrt(.qf.non0(object@V, x)): NaNs produced

## Warning in sqrt(.qf.non0(object@V, x)): NaNs produced
pp <- fitted(model)

newdat <- data.frame(Reef.ID = levels(merged$Reef.ID),
                   newDepth = mean(merged$newDepth, na.rm=T))
pred <- predict(model, newdata = newdat, "probs")
new <- data.frame(Reef.ID=as.character(newdat[,1]), pred)
XY4 <- merge(XY3, new, by="Reef.ID", all=T)

PlotOnStaticMap(KBMap, XY4$Latitude, XY4$Longitude)

XY4 <- merge(XY3, new, by="Reef.ID", all=T)
XY4[XY4 == 0.000] <- 0.0000000001
rownames(XY4) <- XY4$Reef.ID
XY4 <- XY4[, -1]
XY4 <- na.omit(XY4)
apply(XY4, MARGIN=1, FUN=function(reef) {
  floating.pie(xpos = reef["X"], ypos = reef["Y"], 
               x=c(reef["Brown.C.x"], reef["Orange.C.x"], reef["Brown.D.x"], 
                   reef["Orange.D.x"]), radius = 7, 
   col = c("turquoise4", "steelblue1", "tomato", "coral"))
})

legend("bottomleft", legend=c("Brown-C", "Orange-C", "Brown-D", "Orange-D"), fill = c("turquoise4", "steelblue1", "tomato", "coral"))

par(new=T, mar=c(15,22,0,0))
HI <- readOGR("coast_n83.shp", "coast_n83") 
HI <- spTransform(HI, CRS("+proj=longlat +datum=NAD83")) 
plot(HI, xlim=c(-158.3, -157.6), ylim=c(21.35, 21.6), lwd=0.4, col="gray", bg="white")
rect(-157.9, 21.41, -157.75, 21.53)
box()

Raw proportion values for Color Morph and Dominant Symbiont at each reef not discounting the effect of depth.

Mantel Test for Spatial Autocorrelation
reef.dists <- dist(cbind(XY4$Longitude, XY4$Latitude))
dom.dists <- bcdist(cbind(XY4$Brown.C.x, XY4$Orange.C.x, XY4$Brown.D.x, XY4$Orange.D.x))
set.seed(12456)
mantel(dom.dists~reef.dists)
##    mantelr      pval1      pval2      pval3  llim.2.5% ulim.97.5% 
##  0.1916069  0.0210000  0.9800000  0.0270000  0.1247632  0.2573410
PlotOnStaticMap(KBMap, XY4$Latitude, XY4$Longitude)

XY4 <- merge(XY3, new, by="Reef.ID", all=T)
XY4[XY4 == 0.000] <- 0.0000000001
rownames(XY4) <- XY4$Reef.ID
XY4 <- XY4[, -1]
XY4 <- na.omit(XY4)
apply(XY4, MARGIN=1, FUN=function(reef) {
  floating.pie(xpos = reef["X"], ypos = reef["Y"], 
               x=c(reef["Brown.C.y"], reef["Orange.C.y"], reef["Brown.D.y"], 
                   reef["Orange.D.y"]), radius = 7, 
               col = c("turquoise4", "steelblue1", "tomato", "coral"))
})

legend("bottomleft", legend=c("Brown-C", "Orange-C", "Brown-D", "Orange-D"), fill = c("turquoise4", "steelblue1", "tomato", "coral"))

par(new=T, mar=c(15,22,0,0))
HI <- readOGR("coast_n83.shp", "coast_n83") 
HI <- spTransform(HI, CRS("+proj=longlat +datum=NAD83")) 
plot(HI, xlim=c(-158.3, -157.6), ylim=c(21.35, 21.6), lwd=0.4, col="gray", bg="white")
rect(-157.9, 21.41, -157.75, 21.53)
box()

Adjusted values of Dominant Symbiont and Color Morph adjusted for the influence of depth.

Mantel Test for Spatial Autocorrelation
reef.dists <- dist(cbind(XY4$Longitude, XY4$Latitude))
dom.dists <- bcdist(cbind(XY4$Brown.C.y, XY4$Orange.C.y, XY4$Brown.D.y, XY4$Orange.D.y))
set.seed(12456)
mantel(dom.dists~reef.dists)
##    mantelr      pval1      pval2      pval3  llim.2.5% ulim.97.5% 
## 0.11488056 0.08300000 0.91800000 0.16500000 0.04445094 0.19470011
Color and Dominant Symbiont per Bay Area

When considering the effect of bay area on the interaction of color morph and dominant symbiont, the proportion of b colonies dominated by D increases as you move from the north to the south of the bay. The proportion increases almost 6x, yet the small number of colonies makes this non-compelling.

Type=table(merged$ColDom, merged$Bay.Area)
prop.table(Type, margin = 2)
##           
##               Central   Northern   Southern
##   Brown.C  0.48076923 0.43589744 0.42574257
##   Orange.C 0.11057692 0.23589744 0.16501650
##   Brown.D  0.04807692 0.00000000 0.05280528
##   Orange.D 0.36057692 0.32820513 0.35643564
chisq.test(Type)
## 
##  Pearson's Chi-squared test
## 
## data:  Type
## X-squared = 20.668, df = 6, p-value = 0.002104

Color Morph per Reef Type

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  results
## X-squared = 0.08426, df = 1, p-value = 0.7716

Symbiont Community Composition per Reef Type

## 
##  Pearson's Chi-squared test
## 
## data:  results
## X-squared = 2.8371, df = 3, p-value = 0.4174

Symbiont Community Composition per Reef ID

## 
##  Pearson's Chi-squared test
## 
## data:  results
## X-squared = 87.127, df = 72, p-value = 0.1081

Symbiont Dominance and Color Morph Across Depths

Dominant Symbiont by Depth

merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
Dom1 <- subset(merged, !is.na(newDepth) & !is.na(Dominant))
results=glm(Dominant~newDepth, family = "binomial", data = merged)
anova(results, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Dominant
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                       705     942.15              
## newDepth  1   129.01       704     813.13 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
logi.hist.plot(Dom1$newDepth, Dom1$Dominant, boxp = FALSE, type = "hist", col="gray", xlabel = "Depth (m)", ylabel = "", ylabel2 = "")
mtext(side = 4, text = "Frequency", line = 3, cex=1)
mtext(side = 4, text = "C                             D", line = 2, cex = 0.75)
mtext(side = 2, text = "Probability of clade C Symbiont", line = 3, cex = 1)

merged$DepthInt <- cut(merged$newDepth, breaks = 0:13)
merged$Dominant2 <- ifelse(merged$Dom=="C", 1, 0)
results=table(merged$Dominant2, merged$DepthInt)
results
##    
##     (0,1] (1,2] (2,3] (3,4] (4,5] (5,6] (6,7] (7,8] (8,9] (9,10] (10,11]
##   0   149    33    27    23    14     4     3     2     2      1       0
##   1    73    52    83    92    46    27    15    12    16      8       5
##    
##     (11,12] (12,13]
##   0       1       0
##   1       0       1
props <- prop.table(results, margin = 2)
par(mar=c(4, 4, 2, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("red", 0.25), alpha("blue", 0.25)), 
        xlab = "", ylab = "",
        space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("C", "D"), fill=c(alpha("blue", 0.25), alpha("red", 0.25)), inset = c(-.2, 0), xpd = NA)
par(new = T)
par(mar=c(4.2, 4, 2, 6))
results=glm(Dominant~newDepth, family = "binomial", data = merged)
fitted <- predict(results, newdata = list(newDepth=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="Depth (m)", ylab="Probability of D-Dominance")

Bars indicate relative frequency of occurence of clade C vs. D dominance at 1m depth intervals when pooling all samples collected. Overlaid line indicates probability of a colony being D-dominated across depths.

Mixture by Depth

merged$Mixture <- ifelse(!merged$Mix=="C" & !merged$Mix=="D", 1, 0)
merged$Mixture2 <- ifelse(!merged$Mix=="C" & !merged$Mix=="D", 0, 1)
results=glm(Mixture~newDepth, family = "binomial", data = merged)
anova(results, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Mixture
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                       705     973.27              
## newDepth  1   93.723       704     879.55 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
results=table(merged$Mixture2, merged$DepthInt)
results
##    
##     (0,1] (1,2] (2,3] (3,4] (4,5] (5,6] (6,7] (7,8] (8,9] (9,10] (10,11]
##   0   154    44    39    28    18     8     4     2     5      2       1
##   1    68    41    71    87    42    23    14    12    13      7       4
##    
##     (11,12] (12,13]
##   0       1       0
##   1       0       1
props <- prop.table(results, margin = 2)
par(mar=c(4, 4, 2, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("red", 0.25), alpha("blue", 0.25)), 
        xlab = "", ylab = "",
        space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("Mix", "Non-Mix"), fill=c(alpha("blue", 0.25), alpha("red", 0.25)), inset = c(-.23, 0), xpd = NA)
par(new = T)
par(mar=c(4.2, 4, 2, 6))
results=glm(Mixture~newDepth, family = "binomial", data = merged)
fitted <- predict(results, newdata = list(newDepth=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="Depth (m)", ylab="Mixture Proportion")

Bars indicate relative frequency of occurrence of Mixtures vs. Non-Mixtures at 1m depth intervals when pooling all samples collected. The overlaid line indicates the probability of a colony being a mixture across depths.

Color Morph by Depth

merged$Color <- ifelse(merged$Color.Morph=="Orange", 1, 0)
results=glm(Color~newDepth, family = "binomial", data = merged)
anova(results, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: Color
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                       706     979.08              
## newDepth  1   47.199       705     931.88 6.413e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Color <- subset(merged, !is.na(newDepth) & !is.na(Color))
logi.hist.plot(independ = Color$newDepth, depend = Color$Color, type = "hist", boxp = FALSE, ylabel = "", col="gray", ylabel2 = "", xlabel = "Depth (m)")
mtext(side = 4, text = "Frequency", line = 3, cex=1)
mtext(side = 4, text = "Brown                       Orange", line = 2, cex = 0.75)
mtext(side = 2, text = "Probability of Orange Color Morph", line = 3, cex = 1)

merged$Color <- ifelse(merged$Color.Morph=="Orange", 0, 1)
results=table(merged$Color, merged$DepthInt)
results
##    
##     (0,1] (1,2] (2,3] (3,4] (4,5] (5,6] (6,7] (7,8] (8,9] (9,10] (10,11]
##   0   153    37    51    58    24    11     9     1     6      1       1
##   1    69    48    59    57    36    20     9    13    13      8       4
##    
##     (11,12] (12,13]
##   0       1       0
##   1       0       1
props <- prop.table(results, margin = 2)
par(mar=c(4, 4, 2, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("orange", 0.25), alpha("sienna", 0.25)), 
        xlab = "", ylab = "",
        space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("Brown", "Orange"), fill=c(alpha("sienna", 0.25), alpha("orange", 0.25)), inset = c(-.22, 0), xpd = NA)
par(new = T)
par(mar=c(4.2, 4, 2, 6))
merged$Color2 <- ifelse(merged$Color=="0", 1, 0)
results=glm(Color2~newDepth, family = "binomial", data = merged)
fitted <- predict(results, newdata = list(newDepth=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="Depth (m)", ylab="Probability of Orange-Dominance")

Bars indicate relative frequency of b vs. o color morph dominance at 1m depth intervals when pooling all samples collected. The overlaid line indicates the probablity of a colony to be o across depths.

Color Morph and Dominant Symbiont Interaction

Dominant Symbiont per Depth and Bay Area

model=aov(Dominant~newDepth*Bay.Area, data = merged)
Anova(model, type = 2)
## Anova Table (Type II tests)
## 
## Response: Dominant
##                    Sum Sq  Df  F value  Pr(>F)    
## newDepth           25.026   1 125.1913 < 2e-16 ***
## Bay.Area            0.954   2   2.3870 0.09265 .  
## newDepth:Bay.Area   1.555   2   3.8886 0.02092 *  
## Residuals         139.933 700                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Depth at which Orange switches from D to C Dominance

While b was always C-dominated, o was observed to switch from D to C dominance. The depth threshold where this switch occurred is represented here.

threshdepth <- function(color) {
  df <- subset(merged, Color.Morph==color)
  plot(df$Dominant2~df$newDepth, xlab="Depth (m)", ylab = "Probability of Clade C Symbiont",
       main=color)
  abline(h = 0.5, lty=2)
  results=glm(Dominant2~newDepth, family = "binomial", data = df)
  pval <- data.frame(coef(summary(results)))$`Pr...z..`[2]
  mtext(side=3, text=pval)
  newdata <- list(newDepth=seq(0,11,0.01))
  fitted <- predict(results, newdata = newdata, type = "response")
  lines(fitted ~ seq(0,11,0.01))
  thresh <- ifelse(pval < 0.05,
                   newdata$newDepth[which(diff(sign(fitted - 0.5))!=0)], NA)
  return(thresh)
}
sapply(levels(merged$Color.Morph), FUN=threshdepth)
## list()
merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
dat <- aggregate(data.frame(prop=merged$Dominant), by=list(Reef.ID=merged$Reef.ID), 
                 FUN=mean, na.rm=T) #prop D/reef
mod <- glm(Dominant ~ newDepth + Reef.ID, data=merged)
mod2 <- glm(Dominant ~ newDepth * Reef.ID, data=merged)
library(lsmeans)
lsm <- lsmeans(mod, specs="Reef.ID")
lsm2 <- lsmeans(mod2, specs="Reef.ID")
## NOTE: Results may be misleading due to involvement in interactions
res <- merge(dat, data.frame(summary(lsm))[,c(1:2)], by="Reef.ID")
res <- merge(res, data.frame(summary(lsm2))[,c(1:2)], by="Reef.ID")
colnames(res) <- c("Reef.ID", "raw.c", "c.adj", "c.adj.int")

XY2 <- merge(XY, res, by = "Reef.ID")
XY2$d.adj <- 1-XY2$c.adj
XY2$d.adj.int <- 1-XY2$c.adj.int
XY2$raw.d <- 1-XY2$raw.c

manteltable = table(merged$Dom, merged$Color.Morph, merged$Reef.ID)
nc <- aggregate(interaction(merged$Color.Morph, merged$Dom), 
                by=list(merged$Reef.ID), FUN=table)
nc <- data.frame(Reef.ID=as.character(nc[,1]), prop.table(nc[,2], margin=1))
XY3 <- merge(XY2, nc, by="Reef.ID", all=T)

merged$ColDom <- interaction(merged$Color.Morph, merged$Dom)
model <- multinom(ColDom~newDepth+Reef.ID, merged)
## # weights:  108 (78 variable)
## initial  value 978.723819 
## iter  10 value 700.001242
## iter  20 value 684.814917
## iter  30 value 683.516711
## iter  40 value 683.175651
## iter  50 value 683.081045
## iter  60 value 683.060482
## iter  70 value 683.058377
## final  value 683.058329 
## converged
means <- lsmeans(model, specs = "Reef.ID")
means
##  Reef.ID     prob           SE df lower.CL upper.CL
##  10          0.25 4.032745e-10 78     0.25     0.25
##  13          0.25          NaN 78      NaN      NaN
##  14          0.25          NaN 78      NaN      NaN
##  18          0.25          NaN 78      NaN      NaN
##  19          0.25 1.425791e-10 78     0.25     0.25
##  20          0.25          NaN 78      NaN      NaN
##  21          0.25          NaN 78      NaN      NaN
##  25          0.25          NaN 78      NaN      NaN
##  26          0.25 1.840688e-10 78     0.25     0.25
##  30          0.25          NaN 78      NaN      NaN
##  38          0.25 2.328306e-10 78     0.25     0.25
##  42          0.25          NaN 78      NaN      NaN
##  46          0.25          NaN 78      NaN      NaN
##  5           0.25 4.136431e-10 78     0.25     0.25
##  Deep        0.25          NaN 78      NaN      NaN
##  F1-46       0.25 2.851581e-10 78     0.25     0.25
##  F2-25       0.25          NaN 78      NaN      NaN
##  F3-18       0.25 2.328306e-10 78     0.25     0.25
##  F4-34       0.25 2.328306e-10 78     0.25     0.25
##  F5-34       0.25 3.292723e-10 78     0.25     0.25
##  F6-Lilipuna 0.25 4.197414e-10 78     0.25     0.25
##  F7-18       0.25 2.469542e-10 78     0.25     0.25
##  F8-10       0.25          NaN 78      NaN      NaN
##  F9-5        0.25 1.164153e-10 78     0.25     0.25
##  HIMB        0.25          NaN 78      NaN      NaN
## 
## Results are averaged over the levels of: ColDom 
## Confidence level used: 0.95
pp <- fitted(model)

depthmod <- multinom(ColDom ~ newDepth, merged)
## # weights:  12 (6 variable)
## initial  value 978.723819 
## iter  10 value 745.177967
## final  value 745.142241 
## converged
depthdat <- data.frame(newDepth = seq(0,12,0.5))
depthpred <- predict(depthmod, newdata=depthdat, "probs")
#plot(NA, xlim=c(0,12), ylim=c(0,1))
stackpoly(depthpred[,c(4,2,3,1)], stack=T, col=rev(c("burlywood4","burlywood3","darkorange3","darkorange")),
          ylim=c(0,1), xaxlab=depthdat$newDepth, xlab="Depth", ylab="Probability")
legend("top", legend=c("Brown.C",  "Brown.D", "Orange.C","Orange.D"), pch=22,
       pt.bg=c("burlywood4","burlywood3","darkorange3","darkorange"), inset=-0.3, xpd=NA, bty="n")

df <- subset(merged, Color.Morph=="Orange")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
par(mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="orange", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)
df <- subset(merged, Color.Morph=="Brown")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
fitted <- predict(results, newdata = newdata, type = "response")
lines(fitted~seq(0,11,0.01), col="sienna", lwd=3)
mtext(side = 1, text = "Depth (m)", line = 3, cex = 1)
mtext(side = 2, text = "Probability of Clade D Symbiont", line = 3, cex = 1)
legend("topright", legend=c("Brown", "Orange"), fill=c("sienna", "orange"), inset = c(-.22, 0), xpd = NA)

Color Morph by Depth and Reef Type

model2=aov(Color~newDepth*Reef.Type, data = merged)
Anova(model2, type = 2)
## Anova Table (Type II tests)
## 
## Response: Color
##                     Sum Sq  Df F value    Pr(>F)    
## newDepth            11.490   1 49.1052 5.693e-12 ***
## Reef.Type            0.188   1  0.8039    0.3702    
## newDepth:Reef.Type   0.478   1  2.0448    0.1532    
## Residuals          164.489 703                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
df <- subset(merged, Reef.Type=="Patch")
results=glm(Color2~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
par(mar=c(4, 4, 2, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="dodgerblue3", lwd=3, xlab="", ylab="", axisnames=FALSE)
abline(h = 0.5, lty=2)
df <- subset(merged, Reef.Type=="Fringe")
results=glm(Color2~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
fitted <- predict(results, newdata = newdata, type = "response")
lines(fitted~seq(0,11,0.01), col="brown1", lwd=3)
mtext(side = 1, text = "Depth (m)", line = 3, cex = 1)
mtext(side = 2, text = "Probability of Orange Color Morph", line = 3, cex = 1)
legend("topright", legend=c("Patch", "Fringe"), fill=c("dodgerblue3", "brown1"), inset = c(-.22, 0), xpd = NA)

Figures

Symbiont Community Composition per Dominant Symbiont

results=table(Symcap$Mix, Symcap$Dom)
chisq.test(results)
prop.table(results, margin = 2)
par(mar=c(4, 4, 2, 6))
barplot(prop.table(results, margin = 2), col = c(alpha("blue", 0.75), alpha("blue", 0.25), alpha("red", 0.25), alpha("red", 0.75)), xlab = "Dominant Symbiont", ylab = "Symbiont Mixture Proportion")
legend("topright", legend=c("C", "CD", "DC", "D"), fill=c(alpha("blue", 0.75), alpha("blue", 0.25), alpha("red", 0.25), alpha("red", 0.75)), inset = c(-.2, 0), xpd = NA)

When looking at D-only colonies, the Ct values are on the higher end of the spectrum (35 or greater) indicating poor amplification and the potential for C amplification being pushed back in the cycle. This is supported by the fact that 5 of the 9 D-only colonies had C present in 1 qPCR replicate.

Proportion of Clade D when Present

propD <- merged$propD[which(merged$propD > 0 & merged$propD < 1)]

propDHist <- subset(merged, propD > 0 & propD < 1)
propDHist$propD <- cut(propDHist$propD, breaks = 5)
DCol=table(propDHist$Color.Morph, propDHist$propD)
z <- with(subset(merged, propD==0), table(Color.Morph))
o <- with(subset(merged, propD==1), table(Color.Morph))
DCol <- cbind(z, DCol, o)
par(mar=c(4, 4, 2, 0))
bars <- barplot(DCol, xlab = "Clade D Proportion", ylab = "Number of Samples", main = "", col = c(alpha("sienna", 0.55), alpha("orange", 0.55)), axisnames = FALSE, space = c(0,0.3,0,0,0,0,0.3))
axis(side = 1, at = c(1.3, 2.3, 3.3, 4.3, 5.3, 6.3), labels = c(">0%", "20%", "40%", "60%", "80%", "<100%"))
axis(side = 1, at = c(0.5, 7.2), labels = c("All C", "All D"), tick = FALSE)

Depth Influence on Dominant Symbiont and Color Morph

par(mfrow=c(3,1))

merged$DepthInt <- cut(merged$newDepth, breaks = 0:13)
merged$Dominant <- ifelse(merged$Dom=="C", 0, 1)
merged$Dominant2 <- ifelse(merged$Dom=="C", 1, 0)
results=table(merged$Dominant2, merged$DepthInt)
results
props <- prop.table(results, margin = 2)
par(mar=c(2, 4, 2, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("red", 0.25), alpha("blue", 0.25)), 
        xlab = "", ylab = "",
        space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("C", "D"), fill=c(alpha("blue", 0.25), alpha("red", 0.25)), inset = c(0, 0), xpd = NA)
par(new = T)
par(mar=c(2.1, 4, 2, 6))
results=glm(Dominant~newDepth, family = "binomial", data = merged)
fitted <- predict(results, newdata = list(newDepth=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="", ylab="Probability of D-Dominance")

merged$Color <- ifelse(merged$Color.Morph=="Orange", 0, 1)
results=table(merged$Color, merged$DepthInt)
results
props <- prop.table(results, margin = 2)
par(mar=c(3, 4, 1, 6), lwd = 0.25)
barplot(props[,1:11], col = c(alpha("orange", 0.25), alpha("sienna", 0.25)), 
        xlab = "", ylab = "Probability of Orange-Dominance",
        space = 0, xaxs="i", yaxs="i", axisnames = FALSE)
par(lwd=1)
legend("topright", legend=c("Brown", "Orange"), fill=c(alpha("sienna", 0.25), alpha("orange", 0.25)), inset = c(0, 0), xpd = NA)
par(new = T)
par(mar=c(3.1, 4, 1, 6))
merged$Color2 <- ifelse(merged$Color=="0", 1, 0)
results=glm(Color2~newDepth, family = "binomial", data = merged)
fitted <- predict(results, newdata = list(newDepth=seq(0,11,0.1)), type = "response")
plot(fitted~seq(0,11,0.1), xaxs="i", yaxs="i", xlim=c(0,11), ylim=c(0,1), type="l", lwd = 3, xlab="", ylab="")

df <- subset(merged, Color.Morph=="Orange")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
par(mar=c(4, 4, 0, 6))
fitted <- predict(results, newdata = newdata, type = "response")
plot(fitted ~ seq(0,11,0.01), ylim = c(0,1), type="l", col="orange", lwd=3, xlab="Depth (m)", ylab="Probabilty of D-Dominance", axisnames=FALSE, xaxs = "i", yaxs = "i")
abline(h = 0.5, lty=2)
df <- subset(merged, Color.Morph=="Brown")
results=glm(Dominant~newDepth, family = "binomial", data = df)
newdata <- list(newDepth=seq(0,11,0.01))
fitted <- predict(results, newdata = newdata, type = "response")
lines(fitted~seq(0, 11, 0.01), col="sienna", lwd=3)
legend("topright", legend=c("Brown", "Orange"), fill=c("sienna", "orange"), inset = c(0, 0), xpd = NA)

#results=table(merged$Dominant2, merged$Color.Morph)
#chisq.test(results)
#prop.table(results, margin = 2)
#par(new=T, mar=c(10, 10, .5, 6.3))
#barplot(prop.table(results, margin = 2), col = c(alpha("red", 0.35), alpha("blue", 0.35)), xlab = "", ylab = "", yaxt = 'n')
#legend("topright", legend=c("C", "D"), fill=c(alpha("blue", 0.35), alpha("red", 0.35)), inset = c(0, 0), xpd = NA)

Color Morph and Dominant Symbiont at each Reef

par(oma=c(3,3,0,0))
PlotOnStaticMap(KBMap, XY4$Latitude, XY4$Longitude)
axis(1, at = LatLon2XY.centered(KBMap, NA, c(-157.85, -157.81, -157.77))$newX, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("157.85°W", "157.81°W", "157.77°W"), padj = -2.5, cex.axis = 0.75)
axis(2, at = LatLon2XY.centered(KBMap, c(21.42, 21.46, 21.50), NA)$newY, tcl=0.5, line = 0.5, col = "ghostwhite", col.ticks = "black", lwd = 1, outer = TRUE, labels = c("21.42°N", "21.46°N", "21.50°N"), padj = 0.5, hadj = 0.60, las = 1, cex.axis = 0.75)

XY4 <- merge(XY3, new, by="Reef.ID", all=T)
XY4[XY4 == 0.000] <- 0.0000000001
rownames(XY4) <- XY4$Reef.ID
XY4 <- XY4[, -1]
XY4 <- na.omit(XY4)
apply(XY4, MARGIN=1, FUN=function(reef) {
  floating.pie(xpos = reef["X"], ypos = reef["Y"], 
               x=c(reef["Brown.C.y"], reef["Orange.C.y"], reef["Brown.D.y"], 
                   reef["Orange.D.y"]), radius = 7, 
               col = c("royalblue3", "cornflowerblue", "red", "rosybrown2"))
})

legend("bottomleft", legend=c("Brown-C", "Orange-C", "Brown-D", "Orange-D"), fill = c("royalblue3", "cornflowerblue", "red", "rosybrown2"))

par(new=T, mar=c(14,21,0,0))
HI <- readOGR("coast_n83.shp", "coast_n83") 
HI <- spTransform(HI, CRS("+proj=longlat +datum=NAD83")) 
plot(HI, xlim=c(-158.3, -157.6), ylim=c(21.35, 21.6), lwd=0.4, col="gray", bg="white")
rect(-157.87, 21.41, -157.75, 21.52)
box()

Adjusted values of Dominant Symbiont and Color Morph adjusted for the influence of depth.